%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%-----------------------  Simulate Equilibrium  --------------------------%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

%% Setup workspace
close all
clear 

machine = '';

set(0,'defaultTextInterpreter','latex')

M = 15;

%% Load matching function estimates

max_length = 8;     % we drop ships that have been waiting at port for longer than "max_length" weeks
exclude_first = 0;
folder = "/Users/giulia/Dropbox/Matching Ships/coding work/matching function";
path_mf = strcat(folder,'/output/mf_',num2str(max_length),'_',num2str(exclude_first),'.mat');
load(path_mf)

options = optimset('Algorithm','interior-point','display','on','MaxIter',90000,'MaxFunEvals',90000,'TolFun',10^-9,'TolX',10^-9,'TolCon',10^-9);

%
% For the purpose of the counterfactual simulations,
% we approximate the matching function with a Cobb douglas
% 

cobb_par = zeros(15,2);

for i =1:15
    col = freights(:,i)~=0;
    s = ships(col,i);
    f = freights(col,i);
    m = matches(col,i);
    lb = [0,0];
    ub = [Inf,0.95];
    theta0 = [rand(1),rand(1)];
    [cobb_par(i,:),~,~] = fmincon(@(theta)cobb(theta,f,s,m),theta0,[],[],[],[],lb,ub,[],options);
end

scale = 10^5;

s_values = [1,linspace(0.95,1.95,10)];
v_values = linspace(0.75,5,11);

clearvars -except s_values v_values ss hh scale machine cobb_par

%% Load cost estimates

path_ships = '';

path_exports = '';

load(path_ships);

load(path_exports);

g_euler = 0.57721;


%% Data
market_name = {'N America WC','N America EC','Central America',...
'S America WC','S America EC','West Africa','Mediterranean','Baltic States'...
'S Africa','Middle East','India','SE Asia','China','Australia','Japan-Korea'};

%% Simulated Equilibrium with approximate cobb douglas

%
% Initialize the starting point
%
perc = 0.1;

f_0 = f_init;

s_0 = s_init;

U_0 = U_obs;

J_0 = J_obs;

Jf_0 = ones(M,M);

m_0 = zeros(M,1);

%
% Start iteration
%
max_iter = 500000;
iter = 1;
g = 0.45;
G = 3;
eps = exp(-5);
dif = 10;

while dif>eps && iter<max_iter

    %
    % Update the number of matches, given the matching function, and the
    % guesses for exporters f_0 and ships s_0
    %
    for i = 1:M

        gamma_obs = cobb_par(i,2);

        m_0(i) = cobb_par(i,1);
        m_0(i) = m_0(i).*(s_0(i).^(1-gamma_obs));
        m_0(i) = m_0(i).*(f_0(i).^(gamma_obs));

        if m_0(i) >min(s_0(i),f_0(i))
            m_0(i) = min(s_0(i),f_0(i)); 
        end

    end

    lambda_f_0 = (m_0./f_0);

    lambda_0 = (m_0./s_0);

    %
    % Update prices, choice probabilities, and value functions 
    %

    max_iter_inner = 1; 
    
    [U_1,W_1,J_1,meanprices,Gdest,Jf,Gunc] = alg_1(max_iter_inner,cs,cu,c_inv,sigma,M,beta,ksi,g_euler,lambda_0,lambda_f_0,gamma_e,r,delta,U_0,K_obs,zeros(M,1),sigma_exp,zeros(M,M));

    if sum(isnan(U_1))>0

        sprintf('there are nans')

        break

    end

    ccp_sim = ccp_equil(W_1,M);
    
    %
    % Update the number of waiting ships and exporters
    %

    s_1 = zeros(M,1);

    for i = 1:M

        s_1(i) = nansum(ccp_sim(:,i).*(s_0-m_0)) + nansum(Gdest(:,i).*m_0);

    end

    dif_s = (s_1-s_0);

    f_1 = delta.*(f_0-m_0) + E.*(1-Gunc(:,1));

    f_1(f_1<0) = 0;

    %
    % Check convergence
    %

    dif_f = (f_0-f_1)./(1+abs(f_0));

    dif_Jf = (Jf_0-Jf)./(1+abs(Jf_0));

    dif_J = (J_0-J_1)./(1+abs(J_0));

    dif_U = (U_0-U_1);
    
    dif = [dif_s;dif_f;dif_J(:);dif_Jf(:);dif_U(:)];

    dif = max(abs(dif))

    coeff = 1-(1./(G+iter^g));

    f_0 = coeff.*f_0 + (1-coeff).*f_1;

    s_0 = coeff.*s_0 + (1-coeff).*s_1;

    U_0 = coeff.*U_0 + (1-coeff).*U_1;

    J_0 = coeff.*J_0 + (1-coeff).*J_1;
    Jf_0 = Jf;
    iter = iter+1;
end

%
% Save the values
%


U_ineff = U_0;

m_ineff = m_0;

f_ineff = f_0;

s_ineff = s_0;

J_ineff = J_1;

W_ineff = W_1;

Gdest_ineff = Gdest;

Gunc_ineff = Gunc;

meanprices_ineff = meanprices;

Jf_ineff = Jf;

ccp_ineff = ccp_sim;

lambda_ineff = lambda_0;

lambda_f_ineff = lambda_f_0;

ineff = Welfare(Gdest_ineff,Gunc_ineff,ccp_ineff,s_ineff,f_ineff,m_ineff,cu./sigma,cs./sigma,c_inv,K_obs,ksi,r,E,sigma);

